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CROSS-SPECTRAL  ANA!  .IS  OF  ACOUSTIC  SIGNALS 
Allan  L.  Gutjahr,  mathematics  department 
Charles  R.  Holmes,  physics  department 
New  Mexico  Institute  of  Mining  and  Technology 

I . INTRODUCTION 

The  spectrum  of  thunder  contains  a considerable  amount 
of  information  about  the  thunder  signal  as  well  as  information 
about  the  lightning  which  is  presumed  to  be  the  cause  of  the 
thunder.  Holmes,  ^ il-  (1971),  McCrory  (1969),  and  Few 
(1968,  1969),  among  others,  have  analyzed  the  acoustic  signal 
generated  by  thunder.  Although  the  shape  of  the  spectrum  is 
not  completely  knov/n,  various  aspects  of  the  spectrum  are 
discussed  in  the  papers  mentioned  above. 

The  main  purposes  of  this  paper  are  to  discuss  the  use 
of  thunder  for  lightning  location,  to  discuss  estimation  of 
the  quantities  that  occur  in  the  digital  analysis  of  signals 
and  to  present  the  statistical  analysis  of  the  estimated 
quantities. 

Most  of  the  procedures  presented  here  are  available  in 
the  literature,  but  they  are  either  described  rather  generally 
or  in  great  mathematical  detail.  In  either  case  they  are 
often  inaccessible  to  the  practitioner  who  is  not  both  a 
statistician  and  electrical  engineer.  In  addition,  the 
available  discussions  are  scattered  throughout  the  various 
publications,  so  that  the  problems  encountered  by  the  user 
are  accentuated  by  differences  in  notation  and  terminology. 


Hence,  this  is  also  an  attempt  to  give  a consistent  treat- 
ment of  all  of  the  aspects  mentioned  above. 

The  following  is  an  overview  and  table  of  contents  of 
the  paper. 

I.  Introduction. 

II.  Review  of  past  work.  This  section  will  include  a 
discussion  of  the  physical  models  which  form  the 
basis  for  the  work  co'ried  out  in  lightning  loca- 
tion and  spectral  analysis.  In  addition,  a 
critique  of  the  past  work  in  lightning  location 
will  be  included. 

III.  Cross-spectral  analysis  and  lightning  location. 
This  section  will  contain  an  exposition  of  the 
method  of  cross-spectral  analysis  and  its  use  in 
lightning  location.  The  location  procedures  and 
confidence  intervals  will  also  be  discussed. 

IV.  Application  of  cross-spectral  analysis  to  location 
of  prima-cord  and  c-4  shots. 

V.  Application  of  cross -spectra  1 analysis  to 
lightning  location. 

The  appendices  will  include  more  detailed  explanations 
of  the  procedures  of  sections  II  and  III.  In  parti cul ar, they 
will  contain  discussions  of  the  various  computer  routines 
used,  and  further  explanations  of  the  formulas  and  procedures 

The  appendices  are: 

A.  Spectral  Analysis 

B.  Cross-spectral  Analysis 


Variance  and  Covariance  Calculations 
Spherical  Location  Procedure 
Plane  Location  Procedure 

Additional  Remarks  on  Spectra  and  Covariances 
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II.  REVIEW  OF  PAST  WORK  AND  THE  PHYSICAL  MODEL 

The  phenomenon  of  lightning  and  its  relation  to  thunder 
is  discussed  extensively  in  Uman  (1969),  Few  (1968),  and 
McCrory  (1969).  We  will  present  a brief  review  of  the  models 
proposed  by  these  authors,  with  particular  emphasis  on  those 
aspects  that  are  relevant  to  our  work  in  lightning  location. 

Thunder  is  associated  with  the  lightning  path.  Various 
theories  have  been  offered  as  to  how  this  association  occurs. 
One  modern  theory  offers  the  explanation  that  thunder  is  due 
to  rapid  heating  in  the  lightning  channel.  This  rapid  heating 
causes  a shock  wave  which  propagates  radially  outv/ard.  As  the 
shock  wave  moves  out,  the  shock-front  pressure  decreases  and 
eventually  the  wave  becomes  a sound  wave,  which,  after  modifi- 
cation by  the  environment,  becomes  the  audible  thunder. 

The  lightning  channel  from  which  the  thunder  emanates 
is  formed  by  a s tepped-1 eader  process.  The  name  stepped-leader 
is  derived  from  the  fact  that  as  the  lightning  channel  is 
built  from  a cloud  to  the  ground,  say,  it  seems  to  move  in 
steps  of  varying  length,  with  pauses  between  the  steps.  The 
average  step  length  is  about  50  meters  and  the  average  duration 
of  each  step  is  about  50  p-sec.  Usually  the  resulting  channel 
Is  branched.  After  the  channel  reaches  the  ground  a rapidly 
moving  return  stroke  occurs,  moving  from  the  ground  to  the 
cloud  and  propagating  along  the  branches.  There  may  be 
subsequent  return  strokes  associated  with  dart-leaders.  The 
dart-leaders  move  from  the  cloud  to  the  ground  along  the  main 
lightning  channel  and  the  ensuing  return  strokes  also  follow 
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only  the  main  channel.  It  is  the  return  strokes  which  cause 
the  rapid  heating,  and  consequently  the  thunder,  in  the  model 
described  above.  Bec-use  the  return  stroke  moves  very 
rapidly,  the  source  of  the  thunder  occurs  almost  simultaneously 
along  the  entire  channel. 

The  description  above  is  for  the  simplest  discharge  to 
ground  (called  a discrete  stroke  - a discrete  flash  contains 
several  such  discrete  strokes).  There  are  also  strokes 
(known  as  continuous  strokes)  where  there  is  a continuing 
current  to  ground  immediately  after  a return  stroke  and  corres- 
ponding hybrid  flashes,  which  consist  of  series  of  continuous 
strokes . 

Few  (1969)  puts  forth  a "string-of-pearl s"  model  of 
lightning,  wherein  the  lightning  channel  supposedly  behaves 
as  a string  of  cylindrical  sources  of  acoustical  energy.  Each 
cylindrical  source  will  appear  as  a spherical  source  at  large 
distances,  where  the  spherical  source  has  a radius  such  that 
it  emits  the  same  energy  as  the  original  cylindrical  source. 

From  the  theory  for  spherical  shock  waves  one  can  then 
obtain  a relation  between  distance  and  the  dominant  frequency 
(McCrory  (1971)],  and  consequently  one  can  predict  where  the 
peaks  in  the  spectral  density  of  thunder  should  occur. 

Various  studies  of  the  power  spectra  of  thunder  have 
been  carried  out  by  Few  (1968),  Few,  ^ ( 1 967),  McCrory  (1971) 

and  Holmes,  et  a1.  (1971).  The  earlier  studies  of  Few,  ^ al . 
(1967),  indicated  a peak  in  the  spectrum  at  about  200  HZ, 
but  more  accurate  estimates  of  the  spectrum  by  McCrory  (1971) 
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and  Holmes,  ^ (1971)  show  that  the  peaks  are  lower  than 

this  --  for  ground  flashes  they  measured  peaks  (after  correc- 
tion for  wind-noise)  in  the  40  to  80  HZ  range.  Some  attenua- 
tion occurs  due  to  propagation  through  the  medium,  with  the 
higher  frequency  portions  of  thunder  suffering  the  greatest 
attenuation. 

Accounting  for  the  attenuation  and  the  non-stationary 
nature  of  the  spectruui  of  thunder,  (i.e.  its  time  dependence) 
Holmes,  (1971)  and  HcCrory  (1971),  found  that  the  power 

spectrum  had  peaks  at  low  frequencies  and  that  these  peaks 
were  inconsistent  with  the  theory  that  thunder  was  caused 
totally  by  the  acoustic  mechanism  described  above.  As  a 
possible  explanation,  these  investigators  put  forth  the  theory 
that  a second  mechanism  was  in  operation  --  this  mechanism 
was  proposed  originally  by  Wilson  (1920)  and  studied  further 
by  Colgate  (1967)  and  Colgate  and  McKee  (1969).  This  accounts 
for  low  frequencies  by  proposing  an  electrostatic  reaction 
due  to  the  collapse  of  the  region  of  charge  storage  within 
the  cloud  at  the  time  of  lightning  discharge.  This  indeed 
would  account  for  the  low  frequency  content  of  the  spectrum 
and  also  is  in  accord  with  the  fact  that  low  frequencies  were 
observed  in  the  later  part  of  the  thunder  signal. 

The  lightning  channel  is  quite  tortuous,  even  on  a 
small  scale.  This  tortuosity  can  introduce  additional 
ripples  on  the  spectrum  as  proposed  by  McCrory  (1971).  In 
addition,  the  fact  that  the  received  signal  is  the  sum  of 
several  signals  from  different  portions  of  the  channel  can 
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also  introduce  ripples  on  the  spectrum.  McCrory  estimates 
this  latter  filtering  effect,  which  he  calls  phase  noise,  and 
concludes  that  the  tortuosity  of  the  channel  is  more  important 
than  the  time  structure  in  explaining  the  hash  or  ripples 
that  appear  on  the  spectrum. 

The  exact  nature  of  a lightning  channel  is  still  not 
known.  Some  authors  have  claimed  that  the  horizontal  portion 
of  the  channel  which  is  within  the  cloud  can  be  quite  large 
in  comparison  to  the  vertical  portion  that  is  observed  (Uman 
(1969)).  Few  has  attempted  to  reconstruct  the  channel  by 
recording  thunder  signals  at  different  locations  and  then 
tracing  back  to  the  source  by  using  co-variance  analysis  on 
the  signals.  Since  our  objective  is  to  do  a similar  recon- 
struction and  since  the  method  relies  rather  heavily  on  Few's 
work  we  next  give  an  extensive  review  of  his  procedure. 

By  considering  the  hydrodynami cal  conditions  right 
after  a lightning  discharge.  Few  obtained  the  equations 

(•—  + y*V)p  + p(V*v)=0  (1) 

(|y  + V • V)c  + gz  + VP  = 0 (2) 

(|^  + v • V)P  + yP(v  • v)  = 0.  (3) 

p is  the  density,  P the  pressure,  v the  velocity  vector  of 
the  sound  wave,  z is  a unit  vector  in  the  z (vertical)  direc- 
tion, g the  force  of  gravity  and  y is  the  ratio  of  the  specific 
heat  of  the  gas  at  constant  pressure  to  the  specific  heat  at 
constant  volume.  The  equations  (1),  (2),  and  (3)  correspond 
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respecti vely,  to  conservation  of  mass,  momentum  and  energy. 

If  the  quantities  v,  p and  P are  written  as  a steady 
state  term  plus  a perturbed  term  we  have,  after  renormaliz- 
ing , 

v(r,t)  = VQ(r)  + c (r ,t)  (4) 

p(r.t)  = pQ(r)[l  + p^(r.t)]  (5) 

P(r,t)  = pQ(r)[l  + P^(r,t)].  (6) 

If  (4)  through  (6)  are  substituted  into  the  previous 
equations  and  if  the  perturbed  terms  are  set  equal  to  0, 
one  obtains 

^ft  ''O  * ^^*^1  ' (c/H)(v^  • z)  + c V • v^  = 0 (7) 

(|y  + Vq  • V)v^  + (g/c)p^  2 + (c/y)V  • 

- c(yH)‘‘'  = 0 (8) 

(ft  + Vg  * V)P^  - (c/H)v^  ♦ 2 + (c/y)^  • = 0.  (9) 

c is  the  .i:!iabatic  s^ecd  of  sound  and  H is  the  scale  height 
of  the  atmosphere.  These  results  are  obtained  by  first 
solving  the  steady  state  equations. 

Finally  expanding  p.|  , P^  and  v.j  in  wave-number  soace 
we  have 


Pi (r,t) 

OD 

= /// 

• OD 

- at] 

dk^dk^dks 

(10) 

Pi (r,t) 

00 

= /// 

• 00 

- Ojt] 

dk.|dk2dk2 

(n) 

Vi (r.t) 

oo 

= /// 

- U)t] 

dk^dk2dk2. 

(12) 
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In  equations  (10}  - '12)  u is  the  angular  frequency.  If  these 
equations  are  substituted  into  equation  (7)  through  (9),  4 
homogeneous  equations  involving  ^10*  o^Q  and  horizontal  and 
vertical  components  of  v^g  are  obtained.  If  this  system  is 
to  have  a unique  solution,  the  determinant  of  the  coefficients 
must  be  zero.  The  resulting  determinant  involves  co,  k and 
“9  ^ (y  - l)g/(YH).  If  cjg,  which  is  quite  small,  is  ignored, 
one  obtains  the  result  that  k • k.  Finally,  ignoring 

gravity  terms,  one  has 

v^Q  - k ikj  P^g  Y • (13) 


Consequ* ntly , P^g  completely  determines  the  solutions  of 
equations  (7)  through  (9). 

fact  PiQ{k,t),  where  time  dependence  is 
permitted,  could  be  estimated  by  using  an  array  of  microphones 
to  obtain  the  thunder  signals  at  several  locations  simulta- 
neously and  then  using  a 4-dimensional  Fast  Fourier  Transform 


to  calculate  Pig(k,t).  This  mightbe  of  some  use  since  P|Q(k,t) 
could  be  used  to  study  the  dependence  of  the  spectrum  on 
distance. 

In  place  of  (10)  through  (12)  one  could  also  assume  that 
j f and  v-j  are  spatially  and  temporally  statistically 
homogeneous  (Tatarski  (1961))  and  then  obtain  results  similar 
to  those  above  by  using  a Four i er-Sti el t jes  representation. 

In  this  case,  P.|g  would  be  a function  of  k and  t. 

The  additional  assumption  that  the  pressure  wave  1s  a 


plane  wave  yields 
v^Q(r,t)  = k 


lik!! 


{Po('')/Pq(»')  ) P^(r,t). 


(14) 
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Then  wri^-ing  P^j^(k)  in  three  mutually  perpendicular  coordin- 
atcj  where  the  coordinate  system  is  oriented  so  that  a unit 
vector  is  along  kp  the  first  coordinate  of  k,  we  have 

Pio(k)  = 

and 

Pl(r.t)  = / Pp(k^)e’'^‘^'‘‘^''0^  dkp  (16) 


Here  6 denotes  the  delta  function  (that  is,  the  coordinate 
system  is  oriented  in  the  direction  of  the  plane  wave). 

Few  then  assumes  that  the  cross-covariance  between 
P.|(rpS)  and  P.j(r2,t)  only  depends  on  the  differences  t-s 
and  f^-rp  which  is  an  assumption  of  cross - stati onari ty . 

This  assumption  is  not  explicitly  stated  by  Few  but  indeed 
it  is  the  heart  of  the  cross-covariance  technique.  If 
P^(r2,t)  = P.j(rpt  + t),.  (i.e.  if  the  pressure  wave  at 
lags  the  pressure  wave  of  r^  by  t units)  then  one  can  estimate 
T,  which  is  assumed  co  be  constant,  by  calculating  the  cross- 
covariance between  P^Cr^.t)  and  P^(rpt)  and  using  the  point 
where  the  cross-covariance  is  maximum  as  the  desired  estimate 
of  the  time  lag. 

The  planarity  of  the  wave  and  the  geometry  of  the 
situation  lead  to  the  following  equation  for  a,  the  direc- 


tion between  r-  - r 


cos  a = 


2 ■ ' 1 
CT 


^0  ' 


and  k ; 


- --l  - ’'0  '3 


(17) 


Vq  is  the  wind  velocity  and  usually  is  neglected. 

If  this  method  is  applied  to  signals  received  at  three 


microphones  placed  in  a triangular  array,  and  if  lightning 
follows  the  s tri ng-of-pearl s model,  then  the  wave  can  be 
traced  back  to  find  the  location  of  the  source.  An  additional 
measurement  of  the  time  of  onset  of  the  lightning  is  also 
required  --  this  can  be  obtained  optically  or  by  field-change 
measurements.  Several  segments  of  the  signals  are  used  to 
get  several  points  on  the  lightning  channel. 

The  assumptions  that  the  cross-covariance  only  depends 
on  the  lag  between  the  received  signals  and  that  the  received 
pressure  wave  is  locally  stationary  can  be  used  to  give  a 
simpler  formulatior-  of  this  problem  which  doesn’t  require 
the  wave  to  be  a plane  wave. 

Conceptually,  the  model  includes  the  assumption  that  a 
point  source  (one  of  the  pearls  on  the  string)  generates  a 
spherical  wave  which  is  received  as  P{r^,t)  at  one  location 
and  P(r2»t)  at  another  location.  Then  if  the  wind  velocity, 
Vg,  is  ignored,  and  P(r2,t)  = P(r^,t  + i),  once  again  one 
can  estimate  the  time  lag  t and  use  the  time  lag  to  find  the 
source,  using  either  spherical  waves  or  plane  waves.  This 
procedure  is  discussed  in  greater  detail  in  the  next  section. 

There  are  several  problems  connected  with  this  method. 

The  actual  source  is  not  spherical  but  rather  cylindrical. 

In  addition  the  waves  from  different  locations  can  interfere 
with  each  other  leading  to  some  confusion  in  the  estimated 
time  lag.  In  fact,  if  two  points  on  the  channel  are  equidistant 
from  the  microphone  then  the  signals  from  those  two  poirits  uTc 
received  simultaneously.  There  may  also  be  some  problems  viith 


reflections  or  echos  although  one  might  be  able  to  filter 
the  signal  to  remove  such  effects.  Finally,  as  discussed 
previously,  it  does  not  appear  that  the  shock  wave  is  the 
only  component  of  thunder. 

One  difficulty  with  the  cross-covariance  technique  is 
that  the  estimated  time  lags  corresponding  to  peaks  in  the 
covariance  function  are  very  difficult  to  study  from  a 
statistical  point  of  view.  Not  only  are  covariances  difficult 
to  treat  statistical ly,  but  in  this  method  one  is  actually 
interested  in  the  time-1 ag  where  the  covariance  function  is 
maximum  and  this  time-lag  is  difficult  to  treat  statistically. 

Few  (1970)  discusses  the  errors  involved  in  the  cross- 
covariance procedure  but  he  neglects  the  fundamental  point 
that  the  estimate  of  t,  the  time-lag,  is  a statistical 
quantity  with  its  own  statistical  behavior.  He  really  treats 
his  estimate  as  ’f  it  were  the  true  value  of  t and  only 
discusses  errors  due  to  discretization  of  the  time  scale  which 
really  ignores  ^'art  of  the  problem. 

Teer  (1973)  also  claims  to  treat  the  statistical  problem 
but  he  only  finds  the  center  of  the  flash  and  constructs  an 
elliptical  region  which  will  encompass  the  signal.  The 
validity  of  this  procedure  is  questionable  since  once  again 
the  statistical  behavior  of  the  estimates  is  not  accounted  for 

If  we  examine  the  simple  model  where  P{r2,  t + t)  = 
P(r.|,t)  and  where  the  signals  are  stationary,  then  we  can 
estimate  t from  the  cross-spectral  density  of  the  two  signals. 
This  cross-spectral  density  is  the  Fourier  transform  of  the 
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cross-covariance  function.  Then,  t can  be  estimated  by 
examining  the  angle  whose  tangent  is  the  ratio  of  the  real 
and  imaginary  parts  of  the  cross-spectral  density.  As  we 
will  see  in  the  next  section  this  arctangent,  as  a function 
of  the  frequency,  v,  has  the  form  x v,  so  standard  regression 
techniques  can  be  used  to  estimate  x.  In  addition,  in  the 
frequency  domain  it  is  easier  to  study  the  statistical 
properties  than  in  the  time  domain.  Consequently,  one  can 
establish  confidence  limits  for  the  estimate  x,  and,  by  using 
propagation  of  error,  confidence  limits  for  the  source  loca- 
tion. In  addition  the  necessary  filtering  can  often  be 
performed  more  readily  in  the  frequency  domain.  For  these 
reasons  we  have  selected  the  cross-spectral  procedure  and 
this  procedure  will  be  discussed  in  the  next  section  and  in 
the  appendices. 
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III.  CROSS-SPECTRAL  ANALYSIS  AND  LIGHTNING  LOCATION 

The  discussion  in  this  section  assumes  that  the  thunder 
signal  is  recorded  at  three  locations,  where  the  locations 
of  the  recording  microphones  form  a triangular  array,  (see 
Figure  1).  Microphone  i is  located  at  position  = (D  i 1 ’'^i2’ 
D^.^)  and  the  signal  received  at  microphone  i is  X^(t). 

Assuming  that  the  signals  are  stationary,  the  cross- 
covariance function  between  signal  i and  signal  j is  defined 
as : 

(s)  = cov[X^  (t)  ,Xj(t  + s)]  (18. a) 

= E[X^(t)Xj(t  + s)]  - E[X.(t)]  E[Xj(t  + s)].  (18. b) 

Few's  method  proceeds  by  calculating  the  sample  cross-variance 

function  and  then  letting  t (the  time  lag  between  signals 

X^.(t)  and  Xj(t))  be  that  point  where  the  covariance  function 

is  maximum.  The  model  used  assumes  that  Xj{t  + t)  = X^(t), 

as  discussed  in  the  previous  section  of  this  paper. 

If  is  the  cross-spectral  density  of  signals 

i and  j then,  as  we  see  in  Appendix  B,  we  can  es*imate 

■f.j,-(u)  by  taking  the  Fourier  transforms  of  X.(t)  and  X.(t) 

and  then  multiplying  the  transforms.  Specifically,  let 
00 

Fj(u)  = / exp{2Tiuti}  Xj(t)dt.  (19) 

— oo 

Then 

f^ j (u)  = F*(u)  Fj(u) . (20) 

where  the  asterik  denote  3 CGmi  plex  conjugation.  In  addition 
if  Xj(t  + t)  = X^(t),  a simple  change  of  variables  leads  to 
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Fj(u)  = exp{2TTUTi}  F^(u), 
and  consequently 


(21) 


f|j(u)  = exp{2TTUTi}  |f^^(u)|^.  (22) 

We  next  break  up  ^^j('j)  into  its  real  and  imaginary 
parts : 

fij(u)  = [c,-j(u)  - iq.j(u)]/2.  (23) 

Cij(u)  is  called  the  co-spectrum  and  q.j(u)  is  called  the 

quadrature  spectrum.  From  eq.  (22)  it  follows  that  e..(u), 

^ J 

the  phase  spectrum,  is 


0^j(u)  = arctan(-  q^ j (u)/c (u) ) = 2rTU.  (24) 

Hence,  the  phase  spectrum  should  be  a straight  line 

through  the  origin  with  slope  Z-nz  if  signal  j lags  signal  i 

by  T units.  The  phase  can  be  estimated  by  finding  the  value 

h 

of  T that  minimizes  E [‘0ij(u^)  - Ztttu,^]^.  Appendix  B dis- 
cusses  the  cross-spectral  calculations  in  greater  detail. 

The  precision  of  the  estimate  of  x can  be  ascertained  by 
using  standard  statistical  procedures.  The  necessary  calcula- 
tions are  all  given  in  Appendix  C. 

The  cross-spectral  approach  and  the  cross-covariance 
procedure  both  proceed  by  partitioning  the  initial  signals 
into  smaller  subsignals.  After  the  partitioning,  each 
segment  is  analyzed  as  above. 

The  spectral  and  cross-spectral  estimates  are  all 
smoothed  in  order  to  stabilize  the  variance  of  the  estimates. 

In  addition, one  can  also  estimate  the  cross-spectrum  by 
transforming  the  cross-covariance  function  instead  of  the 
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original  signals.  Once  again  the  details  are  given  in 
Appendices  A and  B. 
a . Spherical  location  routine 

Let  S be  the  time  from  the  initiation  of  the  source 
until  its  reception  at  microphone  1 and  let  t.|2  and  be 
the  time  lags  estimated  between  microphones  1 and  2,  and 
1 and  3 respectively. 

Then,  if  r'  = ( >"1 » >^2 ’ *"3 ^ location  of  the  source 

in  a new  system  of  coordinates  where  microphone  1 has 
coordinate  (0,0,0),  microphone  2 is  on  the  x-axis  of  the 
primed  system  and  the  x-y  plane  of  the  primed  system  con- 
tains all  3 microphones,  we  can  obtain  the  equations 

S?  = (r')^  + (rp2  + (rp2  (25. a) 

= {r\  - x')2  + (rp-  + (r^)2  (25. b) 

^3  " ^'”1  ■ * ^*"2  ■ ^3^  '"3^'  (25. c) 

Here  6^  = (cS)^  6^  = [c(S  + = ^^(5  + 

c is  the  speed  ot  sound,  and  xl  , y^I  are  coordinates,  in  the 
primed  system  of  microphone  i. 

If  equation  (25. b)  is  subtracted  from  equation  (25. a) 
and  similarily  if  equation  (25. c)  is  subtracted  from  equation 
(25. a)  one  can  obtain  r j , r^,  and  r^  in  terms  of  S,  T12  and 

■^13' 

Finally,  these  coordinates  can  be  translated  back  to 
the  original  frame  of  reference  by  noting  that  the  primed 
unit  vectors  can  be  expressed  as  linear  combinations  of  the 
unprimed  unit  vectors.  The  variance  and  covariance  matrix 
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of  the  estimated  source  location  can  be  obtained  by  using 
propagation  of  error.  The  details  of  this  procedure  are 
given  in  Appendix  D. 

If  (1-a)  X 100%  confidence  intervals  are  desired  for  the 
source  location,  approximate  limits  can  be  obtained  by  using 
r^.  ± Cvar(r^)]^^^  ^a/6’  ’ ~ where  f^.  is  the  i^^ 

estimated  coordinate  value,  is  the  upper  a/6  point  of 

the  standard  normal  distribution  and  var(r^)  is  the  estimated 
variance  of  r^ . More  precise  estimates  can  also  be  obtained 
by  using  a chi-square  distribution  and  the  covariance  matrix 
of  the  r-values,  but  the  resulting  confidence  region  is  an 
ellipsoid  rather  than  a box  and  is  somewhat  harder  to  visualize. 
The  above  box-like  regions  are  known  as  Bonferroni  confidence 
regions . 


b . Plane  wave  location 

If  the  received  sound  wave  is  a plane  wave  then  one  can 
estimate  the  velocity  vector  of  the  wave  by  using  the  lags, 
T..,  and  then  once  again  trace  back  to  the  source. 

’ w 

Here  the  direction  cosines  of  the  angles  between  the 
velocity  of  the  wave  front  and  the  vectors  connecting  micro- 
phones 2 and  3 to  microphone  1 are  obtained. 

In  particular  if  v is  the  unit  velocity  vector  of  the 
plane  wave,  and  if  a.,  is  the  angle  between  v and  the  vector 

* w 

connecting  microphone  i and  j,  then 


CT. 


cos  a 


JJ.  = 


T J 


PiPj 


P .P  . 
1 J 


)ivll 


i , j = 1,2,3. 


These  2 equations  along  with  |jvlj  = 1 can  be  used  to  find  v. 

(Here  P .P . is  the  vector  connecting  microphone  i and  j.) 

^ J 

The  location  of  the  source  in  this  case  is  cSv.  Once 
again  covariances  for  the  estimates  are  obtained  by  propaga 
tion  of  error,  after  linearizing  the  resulting  expressions 
in  T-J2  and  Appendix  E includes  the  covariance  calcula 

tions  and  a discussion  of  the  plane  location  procedure. 
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IV.  LOCATION  FOR  C-4  AND  PRIHA-CORD  SHOTS 

The  methods  were  tested  by  using  experimental  explosive 
charges  with  point  sources  (C-4  shots)  and  line  sources  (Frima- 
Cord  explosives).  The  explosives  were  suspended  from  balloons 
and  the  approximate  positions  obtained  by  Theodolite  Survey 
data.  The  tests  were  carried  out  near  Langmuir  Laboratory  in 
the  Magdalena  Mountains  of  New  Mexico.  The  co-ordinates  used 
in  this  section  and  the  next  one  are  all  related  to  the  radar 
tower  at  Langmuir  Laboratory. 

Several  different  micropnone  networks  were  used  in  these 
reconstructions.  The  three  networks  used  are  designated  as 
the  Saddle,  West  Knoll,  and  Solar  Tower  networks,  corresponding 
to  their  locations  in  South  Ba1dy  near  the  laboratory.  However, 
the  spacing  between  microphones  is  not  the  same  for  all  recon- 
structions as  is  noted  below. 

A considerable  amount  of  output  is  generated  by  the  pro- 
grams used  to  do  the  reconstruction  and  consequently  only  one 
fairly  complete  set  is  shown  below. 

C-4  Shots,  Summer  1977: 

Two  microphone  networks  were  used  for  these  reconstruct  ions  - 
the  Saddle  and  the  West  Knoll  network.  Each  network  consisted 
of  three  microphone  stations  in  a triangular  array  with  approxi- 
mately thirty  meters  between  the  stations. 

Figures  2-a  through  Z-t  show  the  output  for  C-4  shot  number 
one  for  the  Saddle  network.  Figure  2-a  is  a graph  of  three 
recorded  microphone  signals  for  the  Saddle  network.  Figure  2-t 
summarizes  some  of  the  information  about  the  signals  and  the 


FIGURE  2- A 
C-4  Number  1,  1977 
Saddle  Network 
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SIGNAL  LFNGTHa  0.200  SEC,  FULL  SCALE  DEUIATION«10 . 000  NEWTON/SQUARE  METER 
START  TIMES: 

SIGttlB  — II >45:20.308 
SIG»2B  — 11 :45>20  273 
SIC»3B  — 11  >45  ■•20. 243 
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FIGURE  2-e 


HERTZ  :l:10  ‘ 

DRY  233.  1977  — SPECTRAL  DENSITY  FOR  SIGNAL  »1B 
NU«  33.  BUa  3?  5.  LEHGTHa  0,200  SEC. 

SIGtlB.  STHRTS  11  ^4'5>20,3O8  MST. 


FIGURE  2-f 


FIGURE  2-J 


OAY  233,  1977  — SOURCE  LOCATION  FOR  EUENT  SPECIFIED  B'. 
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estinc+ion  techniques  used.  Figure  2-c  gives  an  approximation 
to  the  source  location  using  only  the  starting  times  of  the 
signals.  Figure  2-d  indicates  the  cross -corre 1 at  ion  functions. 
Note  that  even  though  (from  Figure  2-3)  we  might  guess  that  the 
signals  are  highly  correlated  the  peaks  on  the  cross -correl a ti on 
functions  are  not  very  high.  The  log  of  the  spectral  density 
at  one  of  the  three  microphones  is  shown  in  Figure  2-e,  along 
with  a 95%  confidence  interval  for  the  estimated  spectrum.  The 
flatness  of  the  spectrum  suggests  that  the  disturbance  is 
approximately  white  noise.  Figures  2-f  through  2-h  are 
coherencies  between  the  signals  all  of  which  peak  at  about  200 
Hertz.  These  coherencies  measure  the  correlation  between  the 
various  frequency  components  of  the  signals.  The  phase  functions, 
with  confidence  intervals  and  a regression  line  used  to  estimate 
the  time  lags,  are  shown  in  Figures  2-i  through  2-k.  Figure 
2-1  shows  the  final  output  indicating  the  source  location  and 
confidence  intervals  for  the  spherical  location  routine.  (The 
time  window  line  is  superfluous  - originally  one-half  the  width 
of  the  length  of  record  used  was  added  to  the  times  before 
tracing  and  the  uncertainty  introduced  by  this  was  included  in 
the  variance  calculation.  However,  it  appears  more  accurate 
estimates  are  obtained  by  not  including  this  increment.) 

Variance  estimates  (or  more  precisely  confidence  intervals)  are 
only  shown  for  the  spherical  wave  solution.  The  locations  are 
summarized  in  Table  1 and  the  computed  location  contains  the 
surveyed  location  with  its  confidence  region.  The  corrections 
of  the  program  do  lead  to  considerable  change  in  the  predicted 
locations  from  the  initial  guess  of  Figure  2-c  which  corresponds 
to  a "by  eye"  alignment. 
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TABLE  1 

1977  C-4  TEST  SHOTS 

Co-ordinates  (Meters) 


C-4  NO.  1 

X 

y 

z 

Survey 

-166.7 

660.5 

812.6 

Saddle 

-171.3  ± 

10.4 

655,5  ± 

12.1 

811.7 

+ 

3.1 

West  Knoll 

-173.0  ± 

16.2 

681.9  ± 

18.8 

818.4 

21  .4 

C-4  NO.  2 

Survey 

-213.5 

744.0 

568.0 

Saddle 

-229.5  ± 

5.4 

726.3  ± 

6.0 

574.5 

4- 

3.6 

West  Knoll 

-218.6  ± 

11.5 

726.1  ± 

6.6 

568.8 

+ 

6.42 

35 

Figures  3-a  through  3-g  show  some  of  the  relevant  data 
for  the  same  event  using  the  West  Knoll  network.  Here  the 
confidence  intervals  on  the  spectra  and  phase  are  narrower 
but  because  of  the  greater  distance  from  the  source,  wider  on 
the  location.  While  the  confidence  region  for  the  source  doesn't 
contain  the  survey  values  it  just  barely  misses  doing  so  in  the 
y - co-ordi nate  only. 

The  results  are  presented  in  Table  1 along  with  another 
1977  C~4  test  shot.  The  confidence  regions  for  the  second 
experiment  don't  contain  the  surveyed  location  with  the  largest 
discrepancy  appearing  in  the  y - co-ordi nate . Since  errors 
also  exist  in  the  Theodolite  measurements  this  lack  of  overlap 
was  not  considered  to  be  a serious  problem. 

Prima-Cord  Event  1973: 

Some  calculations  were  also  made  with  earlier  C*4  and 
Prima-Cord  tests.  Only  a Prima-Cord  test  will  be  discussed 
here  since  the  C-4  results  are  similar  to  those  presented  above. 

For  this  reconstruction  again  two  networks  we <^6  used.  How- 
ever, the  Saddle  network  had  microphone  stations  on  a one  hundred 
meter  spacing  rather  than  a thirty  meter  spacing. 

Some  of  the  output  is  shown  in  Figures  4-a  through  4-c  and 
the  co-ordinate  values  are  given  in  Table  2.  Figure  5 shows  two- 
dimensional  views  of  the  reconstruction.  The  reconstructions 
agree  reasonably  well  with  each  other.  There  was  a kink  in 
the  cord  which  appears  in  cne  reconstruction  but  not  the  other 
Indicating  the  effect  of  the  aspect  angl?  since  this  is  indeed 
a real  kink.  The  maximum  cross-correlations  occur  at  the  ends 


FIGURE  3-a 
C-4  Number  1,  1977 
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and  at  the  kink.  The  smaller  spaced  network  has  signals  tha 
do  not  decorrelate  as  fast  as  this  larger  spaced  network. 

The  agreement  with  the  surveyed  location  is  not  very  good  at 
the  bottom  but  the  bottom  survey  was  not  very  accurate  so 
the  discrepancy  isn't  toe  alarming. 
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TABLE  2 

1973  PRIMA-CORD  SHOT 
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a.  Saddle  Network 

Co-ordinates 

(Meters ) 

Maxim 

X 

y 

2 

Cross-cor 

-631.7  ± 1.5 

376.0  ± 1.4 

67.9 

±1.1 

.69 

-588.0  i 3.3 

323.1  ± 4.1 

135.6 

t .9 

.49 

-600.1  ± 6.2 

331.1  + 2.7 

168.2 

+ 1.6 

.39 

-610.fi  ± 3.3 

334.6  ± 3.4 

200.7 

+ .9 

.31 

-624 . 0 + 4.9 

342.9  ± 8.7 

232.4 

+ 1 .4 

.40 

-630.3  r 12.0 

354.4  + 6,3 

265.7 

+ 3.2 

.21 

-640.5  ± 35.0 

357.2  + 7.6 

299.1 

± 9.5 

.16 

-641.0  + .6 

348.6  ± 1.0 

379.6  ; 

± .12 

.90 

b.  West  Knoll  Network 

-659.4  + 4.3 

380.7  ± 2.8 

67.3  i 

: 14 

.85 

-660.9  ± 7.6 

344.1  + 5.8 

130.1  ± 

17 

.82 

-663. 9+8.4 

332.3  ± 6.9 

208.8  i 

13.4 

. 72 

-661 .6  ± 7.2 

327.3  + 3.0 

269.0  ± 

9.5 

.86 

-646.6  ± 3.7 

331 .8  ± 3.0 

305.0  + 

4.7 

.87 

-636.4  ± 2.2 

324.2  ± 1.6 

373.9  ± 

2.4 

. 99 

c.  Double  Theodolite 

Survey 

Top;  -648 

356 

385 

Bottom:  -723 

385 

85 

(Best  guess) 

49 


V.  LOCATION  FOR  SOME  LIGHTNING  EVENTS 

This  section  includes  a few  reconstructions  for  lightning 
events.  The  primary  purpose  is  to  show  some  reconstructions 
and  also  to  discuss  procedures  and  problems  associated  with 
reconstruction  of  lightning  paths  by  the  techniques  discussed 
in  this  report.  Note  that  while  the  individual  points  are 
connected  these  connections  are  to  a large  extent  guesses  and 
are  made  to  give  some  order  to  the  data.  The  more  tortorous 
parts  of  the  channels  may  indeed  look  a lot  different  than  what 
the  figures  below  show.  The  examples  will  be  presented  first 
with  a few  comments.  The  section  closes  with  some  more 
detailed  comments  on  methods  and  procedures. 

Event  1 - Day  233,  1977 

Figures  6-a  through  6-f  show  some  results  from  a lightning 
event  as  recorded  on  the  West  Knoll  network  (30  meter  spacings). 
This  event  was  tracked  with  both  the  West  Knoll  and  Saddle 
(30  meter  spacings)  networks  and  the  results  are  shown  in  Table 
3.  The  two-dimensional  views  of  the  reconstruction  are  shown  in 
Figure  7 and  a possible  channel  is  indicated.  While  it  is  diffi- 
cult to  completely  reconstruct  the  lightning  path,  the  general 
region  of  the  stroke  is  fairly  clear  with  both  networks  agreeing 
f a i rl y well. 

Event  2 - Day  209,  1975 

This  event  was  recorded  on  three  networks;  Saddle  (100 
meter  spacing).  West  Knoll  (30  meter  spacing),  and  Solar  Tower 
(30  meter  spacing). 
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TABLE  3 

CO-ORDINATES  FOR  EVENT  1 (Day  233,  1977) 
Lightning  at  14:23:59.166 
a.  Saddle  Network 


Start  Time  Co-ordinates  (Meters) 

X y 2 


14:24:07.124 

-3245 

+ 

50 

-51 

92 

69 

770 

14:24:07.409 

-3276 

42 

-330 

+ 

56 

256 

295 

14:24:07.703 

-3221 

+ 

96 

-554 

125 

726 

+ 

223 

14:24:08.572 

-3602 

+ 

34 

-272 

+ 

172 

752 

+ 

141 

14:24:08.772 

-3664 

82 

-274 

4- 

205 

786 

275 

14:24:09.626 

-3713 

* 

35 

-974 

4 

57 

981 

89 

West  Knoll  Network 

14:24:06.00 

-3356 

± 

28 

-348 

± 

20 

290 

+ 

187 

14:24:06.147 

-3334 

± 

31 

-380 

± 

21 

570 

+ 

95 

14:24:06.347 

-3319 

+ 

31 

-425 

± 

36 

785 

4 

72 

*Times  indicate  these  two  points  are  close  to  each  other. 
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TABLE  4 

CO-ORDINATES  FOR  EVENT  (Day  209,  1975) 
Lightning  at  12:47:07.550 


Saddle  Network 

Co-ordinates 

(Meters ) 

Maximum 

Cross-correlation 

X 

y 

z 

-1230 

1645 

226 

.70 

-1397 

1598 

993 

.75 

-706 

1900 

2351 

.62 

West  Knoll 

Network 

-1151 

1705 

19 

.28 

-1332 

1721 

385 

n.a. 

-1415 

1801 

1280 

.76 

-1489 

1814 

1 500 

n.a. 

-1084 

1762 

2135 

. 74 

-803 

1775 

2339 

.54 

-688 

1858 

2395 

. 71 

Solar  Tower 

Network 

-1317 

1660 

111 

.30 

-1387 

1622 

1156 

. 78 

-1545 

1679 

1532 

.69 

-1587 

1635 

1 744 

. 34 

-882 

1731 

2438 

.64 

-344 

1926 

2570 

.60 

n.a.  means  not  available. 
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The  reconstruction  for  this  (and  the  following  events) 
was  done  with  an  earlier  version  of  the  computer  program. 

The  earlier  version  added  half  the  width  of  the  signal  to  the 
times  of  arrival,  as  indicated  above,  and  also  had  in  error 
in  the  confidence  interval  calculations.  Consequently  the 
results  in  Table  4 don't  include  any  error  limits.  The 
graphs  for  this  event  are  in  Figure  8.  The  three  networks 
give  results  that  agree  very  well.  This  appears  to  be  a very 
clean  cloud-to-ground  stroke  with  little  evidence  of  branching. 

Event  3 - Day  213,  1976 

Figures  9-a  and  b show  a signal  at  the  Solar  Tower  network 
(30  meter  spacing)  for  a more  complex  and  long  record.  Figure 
9-a  shows  one-second  of  the  data  and  Figure  9-b  shows  .300 
seconds.  One  can  pick  out  corresponding  features  from  the 
three  signals  quite  easi-ly  and  the  cros  s - correl  a t i ons  between 
various  signal  segments  were  quite  high. 

Figure  10  shows  the  two-dimensional  views  with  a possible 
path.  In  addition  the  relative  maxima  of  the  maximum  correlation 
values  are  shown  along  the  signal  path.  Here  it  would  appear 
that  cross-corrections  are  maximum  at  corner  points  or  branch- 
ing points. 

Event  4 - Day  210,  1975 

Table  5 and  Figures  11-a,  b,  and  c show  a reconstruction 
for  what  would  appear  to  be  a c 1 oud- to-c 1 oud  stroke.  Three 
networks  were  used  as  the  figures  and  table  indicate.  Again 
this  was  done  with  the  early  version  and  consequently  no 


FIGURE  9-a 

Event  3^  Day  213^  1976 
Solar  Tower  Network 
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FIGURE  9-b 

EVENT  3>  Day  213,  1976 
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TABLE  5 

CO-ORDINATES  FOR  EVENT  4 (Day  210,  1975) 
Lightning  at  11:29:12.760 
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a . 


b. 


Saddle  Network 

Co-ordi nates 

(Meters ) 

Maxiirtum 

Cross-correlation 

X 

y 

z 

942 

5529 

3146 

.72 

-732 

4180 

2382 

.29 

1620 

2755 

3299 

.81 

2361 

2199 

4069 

.71 

816 

5303 

2773 

.74 

185 

6011 

2736 

.74 

680 

2549 

6435 

.68 

-302 

4287 

6890 

.55 

-4818 

6899 

4154 

.84 

West  Knoll 

Network 

1404 

2247 

2979 

.63 

1734 

2841 

3131 

.72 

2033 

2672 

2934 

.73 

1884 

2067 

4397 

.60 

1222 

3309 

5964 

.54 

-2532 

5537 

52C7 

.59 

-2028 

5267 

5617 

.60 

-641 

5258 

6065 

.32 

Solar  Tower 

Network 

-761 

2361 

2809 

.52 

-453 

3901 

2785 

.55 

1457 

3754 

3185 

.50 

1685 

2513 

4360 

. 75 

911 

3540 

5980 

.67 

-2503 

5350 

5447 

.45 

-5805 

5521 

31  37 

.62 

-5055 

7003 

3933 

.76 

j 

i 

1 


$6 

FIGURE  11-B 

Event  Day  210,  1975 
X ; Saddle 
• : West  Knoll 
o ; Solar  Tower 
(y-  z view) 
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FIGURE  11-c 
Event  Day  210^  1975 
X : Saddle 
• : West  Knoll 
o : Solar  Towers 
(x-  Y view) 
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confidence  intervals  are  given  in  the  table.  The  three  recon- 
structions seem  to  agree  fairly  well  and  once  again  suggested 
patterns  for  the  stroke  are  sketched  in  the  figures.  Again 
higher  correlation  values  seem  to  be  associated  with  turning 
points. 

Comments  on  Liqhtr inq  Channel  Reconstruction 

The  purpose  of  this  report  has  been  primarily  to  review 
and  present  the  techniques  of  lightning  channel  reconstruction 
from  acoustic  data.  The  following  comments  pertain  to  the  use 
of  these  techniques. 

Reconstruction  of  lightning  paths  using  the  acoustic 
technique  is  a rather  time  consuming  process.  It  requires  con- 
stant operator  attention  a..d  is  rather  dependent  on  a good 
initial  guess  regarding  s'gnal  alignment.  A reasonably  clean 
signal  (i.e.,  isolated  from  other  noise  sources)  'is  needed  in 
order  for  the  method  to  succei.a.  The  rc:;ults  presented  here 
and  additional  confirmatory  results  using  independent  radar 
measurements  by  Scymansr.i  (197:^)  do  indicate  that  the  technique 
can,  under  proper  conditions  yield  a reasonably  accurate  re- 
construction of  a channel. 

Proble.'-'S  involving  the  physics  cf  the  process  need  more 
investigation.,  however.  For  ex.^mple  a detailed  study  of  the 
interference  expected  between  the  acoustic  signals  from  differ- 
ent segments  of  the  lightning  path  .should  be  made  to  see  what 
effect  such  interference  has  on  the  reconstruction.  Higher 
correlations  do  appear  to  occur  at  branching  or  turning  points 
suggest  such  interference  or  interaction  is  minimized  at  these 


69 


points.  An  examination  of  the  relation  between  tortuosity  and 
correlation  is  clearly  a subject  that  needs  more  study.  Finally 
the  effect  of  assuming  a spherical  rather  than  a cylindrical 
wave  may  also  introduce  some  bias  into  the  reconstruction. 

The  confidence  intervals  given  in  the  paper  are  predicated 
on  an  accurate  model  and  pertain  primarily  to  estimation  tech- 
niques. The  fact  that  for  some  of  the  test  cases  these  inter- 
vals did  not  contain  the  surveyed  point  does  indicate  that  the 
model  may  be  biased  although  not  to  a very  strong  degree. 
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APPENDIX  A 

Stochastic  Analysis  and  Spectral  Analysis 

The  general  time  series  model  deals  with  second  order 
stationary  processes,  which  we  define  below. 

Definition:  X(t)  is  a second  order  stationary  stochastic 

process , i f 

(a)  E(X(t))  = y 

(b)  C0V(X(t)),X(t  + s))  - r(s). 

Thus,  the  process  must  have  a constant  mean  and  a covariance 
function  that  only  depends  on  the  time  difference.  We  use 
E to  denote  expected  value  and  COV  for  covariance. 

C0V(X(t),X(t  + s))  = £[X(t)  • X(t  + s)]  - E[X(t)]  E[X(t  + s)] 

(A.1) 

in  case  X(t)  is  a real  process. 

If  Z and  W are  complex  random  variables  and  if  W*  is  the 
complex  conjugate  of  W*.  then 

C0V(Z,W)  = E(ZW-)-E{Z)  E(W*).  (A. 2) 

Condition  (a)  in  the  definition  above  is  not  as  stringent 
as  condition  (b).  In  particular,  if  E(X(t))  = y(t),  then 
using  X(t)  - u(t)  in  place  of  X(t)  will  eliminate  the  depen- 
dence on  t.  In  addition  E(X(t)  - u(t))  = 0,  and  consequently 
we  often  assume  without  loss  of  generality  that  E(X(t))  = 0. 

For  a second  order  stationary  process, 

VAR(X(t))  = r{0) 

and,  since  VAR(X(t  + s)  - X(t))  ^ 0,  r(0)  > r(s)  for  all  s. 


In  addition  r(s)  ~ r(-s).  Hence,  cos(s)  Is  a valid  covariance 
function  but  s1n(s)  Is  not. 

If  a process  Is  : ot  second-order  stationary,  very  little 
can  be  done.  (See  Ha.indn,  p.  77).  Most  of  the  theory  for 
second-order  stationary  processes  relies  on  the  following 
two  representation  theorems. 

THEOREM  I.  If  X(t)  Is  a second-order  stationary  process  then 
there  exists  a complex  random  process,  Z('.)),  such  that 

X(t)  = / e Z(dv)  (A. 3) 

» 00 

where  COV(Z(Avi),  = 0 if  and  Av^  are  two  disjoint 

Intervals  and  Z(-Av)  = Z*(Av).  (Z(A)  = Z{b)  - Z(a)  if  A = 

[b,a]  is  an  interval  on  the  v axis). 

THEOREM  II.  If  r(s)  is  ttie  covariance  function  of  the  second- 
order  stationary  process  X(t),  then 

00 

r(s)  = / e F(dv)  (A. 4) 

« 00 

where  F(Av)  > 0 and  E(!Z(Av)|^)  = F{Av). 

A proof  of  these  theroems  can  be  found  in  Rosenblatt, 
Chapter  VII  (1963). 

F(v)  is  called  the  spectral  distribution  function.  If 
F(v)  is  differentiable  then  f(v)  = ^ (v)  is  called  the  spectral 
density  for  the  process. 

The  v's  correspond  to  frequencies  while  the  Z(v)'s  assign 
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complex  weights  to  frequency  v.  By  Theorem  I,  a second-order 
stationary  process  is  a mixture  of  sines  and  cosines  with 
random  amplitudes  at  each  frequency.  The  condition  that 
Z(-Av)  = Z*(Av)  ensures  a real  value  for  X(t).  The  usefulness 
of  Theorem  I is  due  to  the  uncorrelated  nature  of  the  random 
amplitudes,  Z(Au). 

In  Theorem  II,  the  definition  of  F(Av)  and  the  fact  that 
Z(-Av)  = Z*(A'j)  imply  F(Av)  = F(-Av)  and  hence  f(u)  = f(-v) 
if  the  spectral  density  exists.  We  will  assume  f(v)  always 
exists  in  the  rest  of  this  paper. 

From  the  relationship  between  a Fourier  transform  and 
its  inverse,  we  can  obtain  Z(Av)  and  f(v)  by  taking  the 
appropriate  transform  of  X(t)  and  r(s).  In  order  to  make 
these  transform  pairs  more  symmetric  in  appearance  and  also 
to  make  the  results  consistent  with  conventional  engineering 
approaches  and  computer  programs,  let  v = 2ttu.  Then 

09 

Z(Au)  = / exp{-  2:rAuti}  X{t)dt  (A. 5) 

“00 
and 

f(u)  = / r(s)ds.  (A. 6) 

— 00 

(In  practica'  calculation,  only  positive  frequencies  are 
used  and  the  resulting  estimate  is  doubled.) 

The  results  of  Theorem  I and  II  show  that  one  can  obtain 
the  spectral  density  in  2 different  ways.  One  way,  via 
equation  A. 6,  is  to  take  the  Fourier  transform  of  the  covar- 
iance function.  The  other  way  Is  to  use  the  fact  that 
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f(u)  Au  = F(au)  = E(|Z(au)|‘^)  and  hence  use  equation  A. 5 to 
obtain  Z{Au)  and  then  square  the  absolute  value.  Both 
procedures  have  been  used  in  theory  and  in  practice  and  each 
method  leads  to  insight  about  the  behavior  of  spectral 
estimates . 

2 

In  particular,  we  will  use  f{u)Au  = E(|Z(Au)|  ) to  see 
why  the  chi-square  distribution  is  used  to  set  confidence 
intervals  on  estimated  spectra  and  also  to  see  why  smoothing 
is  desirable  when  one  estimates  spectra. 

We  will  start  by  assuming  that  Z(u)  is  a complex 
Gaussian  process.  This  means  {Re(Z(Au^)),  Im(Z(Au^.))) 
i = l,...n  are  jointly  normal.  This  also  implies  X(t)  is 
a Gaussian  process  and,  since  E{X(t))  = 0,  E(Z(u))  = 0.  It 
can  be  shown  (Breiman,  Chapter  9,  or  Hannan,  Chapter  2)  that 


C0V(Re(Z(Au)) , Im{Z{Au)))  = 0 
VAR(Re{Z(Au)))  = 

and  VAR(Im(Z{Au}))  = ^4^^-  . 


(A. 7. a) 
(A.7.b) 

(A.7.C) 


Consequently  the  real  and  imaginary  parts  of  Z{Au)  are 
statistically  independent,  normal  random  variables,  each  with 
mean  0 and  variance  ^ ^4^-^  . Hence 


(A. 8) 


has  a chi-squared  distribution  with  2 degrees  of  freedom. 
In  addition,  if  Au  and  Au ' are  non-overlapping  intervals 


then 
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+ 2|X(Au')r 
F(Au') 


(A. 9) 


will  have  a chi-squared  distribution  with  4 degrees  of 
freedom  since  the  2 terms  correspond  to  independent  chi- 
squared  variables,  each  with  2 degrees  of  freedom. 

If  Z(u)  is  not  a Gaussian  process  one  can  still  obtain 
an  approximate  chi-squared  distribution  as  the  argument 
below  shov/s. 

Suppose  Au  = [a,b)  is  an  interval  in  frequency  space. 

Break  up  the  interval  [a,b)  into  n non-overlapping  intervals 
n-1 

[a,b)  = U [a  + i ( b - a ) , a + ( i -t- 1 ) ( b - a ) v 
i = 0 n n ' 

n-1 

= U A. 
i*0 

n-1 

2(Au)  = I Z(A.).  (A. 10) 

i=0 


If  Z{A.)  = Z(a  + (i+1)  b-a ) - Z(i  + i ( b - a ) and  if 

n n 

F(Ag)  = F(A^)  = FCA^)  ...  = F(A^),  (A. 11) 

2 

where  F(A^)  ==  E(|Z(A^)|  ),  then  if  the  central  limit  theorem 
applies,  Z(Au)  is  approximately  Gaussian  and  the  previous 
results  apply  with  F(Au)  = nF(AQ). 

This  approach  also  shows  why  the  variance  of  the  spectrum 
doesn't  decrease  with  increasing  sample  size.  For  example, 
if  2|Z(Au)|  /F(Au)  has  a X“  distribution  with  2 degrees  of 
freed  0«n  I 
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VAR[2|2(Au)1^/F(au)]  = 4 and 
VAR(|Z(Au)1^)  = F^(Au), 

independent  of  the  sample  size.  If  the  spectrum  is  reasonably 
flat,  one  can  improve  estimates  for  the  spectral  density  by 
averaging  over  several  frequency  bands  and  hence  improve  the 
accuracy  of  the  estimates  without  introducing  any  appreciable 
bias . 

If  the  covariance  function  is  used  ’n  the  spectral 
estimation,  a weighting  function  [called  a lag  window)  l/(t) 
is  used  to  smooth  the  spectra,!  estimate.  N(t)  is  a function 
such  that; 

W (o)  = 1 (A . 1 2 . a ) 

WCt)  = W{-t)  (A.12.b) 

andW(t)=Oifltl>L.  {A.12.c) 

The  smoothed  spectrum  is 

fw('j)  = / W(t)  rCt)e”’^''‘^^  dt.  (A. 13) 

-00 

If  w(u)  is  the  Fourier  transform  of  W(t)  and  f(u)  is  the 
Fourier  transform  of  r(t)  then  by  the  convolution  property 
of  transforms, 

09 

Fy{u)=/  w(v)f(u-y)dv.  (A. 14) 

— 00 

w(u)  is  called  the  spectral  window.  Thus,  one  can  smooth 
either  before  or  after  transformation  of  the  covariance 
function . 

A lag  or  spectral  window  can  be  used  to  reduce 
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the  variance  of  the  estimate  but  it  can  also  lead  to  biased 
and  correlated  estimates.  (See  Jenkins  and  Watts,  p.  247). 

If  £ u _<  U2. 

,2  / -i  “ 

C0V(fy(u^)  f„(ii2))  = / W{u^-v)  {W(u2+v)  + W(u2-v)}dv 

* OD 

.....  (A. 15) 


so  the  2 esti mates  are  not  independent.  Neolecting 

oa 

/ W(v)  W(v  + 2v)dv  which  is  small  compared  to  the  integral 

-00  t \ ® 

of  the  square,  VAR(fj^(u))  = - ■ ' / W^(v)dv. 

— 00 

2T  is  the  length  of  time  over  which  the  signal  is 
recorded.  Let 


T 

I = / 


W^{t)dt. 


(A.  16) 


-T 

2 

Then,  VAR(fj^(u))  = f (u)  I/T  and  so  I/T  is  the  proportional 
reduction  in  the  variance  due  to  smoothing  when  the  Ian  window 
W(t)  is  used. 

A concept  occasionally  used,  is  that  of  band-width  for  a 
spectral  window.  If 

w(u)  = 1/b,  -b/2  1 u < b/2 

is  a rectangular  window,  it  has  a unique  band-width,  b,  in  the 
frequency  domain.  The  variance  of  an  estimator  that  is  smoothed 
via  this  window  is  f^(u)/Tb. 

For  non-rectangular  windows,  the  equivalent  band-width 
is  defined  as  b = l/I  and  the  standardized  band-width  is  bT. 

Note  that  the  variance  times  band-width  is  constant  and  hence 
large  band-width  implies  small  bias.  Some  common  windows  and 

and  their  properties  are  discussed  in  Jenkins  and  Hatts, 

Chapter  6. 
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A commonly  used  window  is  the  Tukey  window; 

W(t)  = -1(1  + cos(Tit/L))  lt|  < L.  (A. 17) 

This  lag  window  has  the  associated  spectral  window 

w(u)  = [L  sin{2;ruL)3/(2T:uL)(l  - (2uL)^).  (A. 18) 

For  the  Tukey  window, 

I = / w^(u)  du  = 3L/4  (A. 19. a) 

-00 

d.f.  = 2T/I  = 8T/3L  (A.19.b) 

band-width  = 4/3L  (A.19.c) 

varianceratio=3L/4T.  (A.19.d) 


For  the  truncated  estimate  the  Tukey  spectral  window 
corresponds  to  weighted  averages,  with  weighting  factors  of 
(1/4,  1/2,  1/4)  for  u - Au,  u and  u + Au  respectively. 

(A-i)  Direct  Spectral  Estimates 

Since  most  observations  are  digitized  the  integrals  are 
replaced  by  sums  in  the  calculations.  Suppose  X(t)  is 
observed  at  discrete  times  kAt,  k = -ni,  - m+  l,  ,...m-l 
(n  = 2m  observations),  then 


(m-l)At  . 

/ X(t)e’2’^“^  dt  * 

-mAt 


m-1 

I X{kAt)  {cos  (2::ukAt ) + i s i n ( 2TTkuAt ) } At. 

k = -iii 


If  u,^  = k/nAt,  Au  = 1/nAt,  k = -m, 


...  m , and  if 
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m-1 

~ ^ X(kAt)  cos(2TTukAt)  (A. 20. a) 

k = -tn 

m-1 

B(u)  = I X(kAt)  s i n ( 2TTukAt ) , (A.20.b) 

k = -m 

then  the  sample  spectrum  at  is 

f(U|^)  = 2[A^(u^)  + B^(u^)]  At/n  (A. 21. a) 

atk=±l,±2,  ...,±  (m-1 ) . 

At  the  end  points,  Uj^  = 0 and  Uj^  = l/(2At), 

f(Uj^)  = CA^{uj,)  + B^(ujj)]  At/n  (A. 21.5) 


with  a corresponding  loss  of  1 degree  of  freedom, 

(u^y  has  a chi-squared  distribution  with  2 degrees  of 
freedom,  for  k = ± 1,  i 2,  ...,  ± (m-1)  and  1 degree  of  freedom 
for  k = 0 and  k = m,  for  a total  of  2(m-l  ) + 2 = 2m  = n degrees 
of  freedom.  (iJote  that,  by  the  symmetry  of  the  spectrum, 
there  are  only  (m-1)  + 1 independent  estimates  rather  than 
2(m)  independent  estimates).  These  n degrees  of  freedom  are 
distributed  over  the  various  frequency  bands.  If  the  true 
spectrum  is  reasonably  constant  over  i adjacent  bands  (say, 
over  u ± i A u,  i = 1,  ...,  £/2)  and  if 

. 

f(u)  = z [f(u  + iAu)  + f(u  - iAu)]/i 
i =1 

i i_  £f  ( u ) . 

AtTtuT  ^ chi-squared  distribution  with  2i  degrees  of 
freedom . 


VAR[f(u)]  = 4i/i^ 


= 4/Z  is  now  decreased.  If  successive 
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bands  are  averaged  without  overlapping,  nfZl  independent 
estimates  are  obtained,  each  with  2t  degrees  of  freedom.  If 
overlapping  occurs,  the  estimates  will  not  be  statistically 
independent. 

In  this  direct  procedure  (known  as  the  Cooley-Tukey 
procedure)  a "fader"  or  "taper"  is  usually  applied  to  the 
original  series  before  transformation  in  order  to  reduce 
bias.  This  is  a function  a(t)  such  that  a(t)  increases 


from  0 to  1 in  [0,S^),  equals  1 on  [SpS^)  and  decreases 


to  0 on  [S^.T).  Usually  a(t  + T/2)  = a(T/2  - t). 

A typical  fader,  with  S-j  = T/10,  is 

'^[1  - cosCrt/S^)]  = sin^(-rt/2  S^,  0 < t < 


a(t) 


1 


Si  1 t < T - 


^^[1  - cos(-T(T-t)/S^  )]  = sin^(-(T-t)/2  S^),  T - < t < 1 

(A. 22) 


The  signal  X(t)  is  first  multiplied  by  a(t},  where  no.v 

we  assume  X(t)  is  observed  on  [0,T].  After  multiplication, 

a(t)  X(t)  is  transformed  and  then  averaged  as  above.  To  see 

the  effect  of  the  "taper",  let  u(x)  = a(t/T)  so  that  u(x) 

is  defined  on  [0,1].  If  Z(u)  is  the  transform  of  a(t)  X(t) 

2 

then  the  direct  estimate  of  the  spectrum  is  ;Z(u)|  (2J^t/n) 

at  frequency  u (where  n.it  = T). 

However , 


EdZlull^)  = 


n 

r T 

k = l 


kll 


f (u) 


rr  ffiji. 

0 


{ A , ? 3 ) 
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Hence,  the  true  spectral  estimate  becomes 

2Atl7(u) l^/(n  u^(x)dr).  (A. 24) 

0 

For  the  cosine  taper  above, 

1 j -la 

/ u'^(x)dx  = 2 / s in^{Tx/ . 2)dx  + .8 
0 0 

= 2(|)(J){4)  + *6  = -875. 

The  smoothed  (averaged)  estimates  are  weighted  chi-squared 
variables,  after  division  by  f(u).  The  distribution  of  these 
estimates  is  approximated  by  a constant  times  a chi-squared 
variable,  where  the  constant  and  degrees  of  freedom  are 
obtained  by  equating  means  and  variances. 

d = 2f^(u)/VAR(f(u))  (A. 25. a) 

c = E(f(u))/d  (A.25.b) 

If  m + 1 final  estimates  are  obtained  the  degrees  of 
freedom  become 

’I 

= 2n  [/  u^(x)dx]^/{m+l ) [/  u^{x)dx].  (A. 26. a) 

0 0 

For  the  cosine  taper, 

1 * -I  « 

/ u^(x)dx  = .8  + 2/  sin“(^)d.x  = .8055  (A.26.b) 

0 0 

The  band  width  is  d/2N. 

The  above  is  true  for  non-overlapping  bands.  If  averaging 
occurs  in  overlapping  bands,  the  degrees  of  freedom  are 

(2L  + 1)[/^  u^(x)dx]^/[/^  u^(x)dx] 

0 0 


{A.26.C) 
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where  2L  + 1 is  the  number  of  elementary  bands  averaged. 

In  either  case,  the  degrees  of  freedom  at  each  end  point  are 

half  the  figures  g’ven  above. 

The  computer  routine  use";  a fast  lourier  transform 

coded  in  a subroutine  called  NLOGN.  If  2(1),  1=1,  N 

is  a (complex)  input  array,  NLOGN,  produces  an  output  array 
N 

Z\K)  = E expiZiri  (K-1  ) ( J-1  ) ] Z(J)  (A. 27) 

J = 1 

IC  = 1,  ...,  N.  The  number  of  entries  in  the  input  and 
output  array  are  equal. 

X(k)  = X(kAt>,  k = 0,1,  ...,  r-1  is  the  original 
data  array,  k = 0,  ...,  n-l  and  N = n. 

The  steps  in  the  comput.ir  analysis  are 

1.  Taper  X(k)  with  the  cosine  taper  by  multiplying  the 
taper  and  the  series.  Load  the  result  into  the  Z array. 

2.  Fill  out  the  Z array  with  zeroes  in  order  to  get 

NN  = 2M(2L  + 1)  - 1 points  if  M + 1 estimates  are  desired. 

3.  Call  NLOGN  tc  get  Z(I),  I-l , 2,  ...,  MN. 

Tne  elementary  spectral  estimaces  are 

F(I)  = 2tZ(i  ) <.t/.875  NN  (A.  28) 

4.  Average  the  el'.-mentary  estimates, 

L 

(a)  2 F(I)/L  is  the  estimate  at  frequency  0. 

i = ■ 

(2K+1 )L+K 

(t)  { E F(I)}/2Li  1 is  the  estimate  at  frequency 

I=(2L-1 )L+K 

(2L  + 1 )K/2At. 

M(2L^'. 

(c)  t F(I)  is  the  estimate  at  frequency  l/2At. 

I = ( 2ii  - 1 ) L+M 


More  closely  spaced  estimates  could  be  obtained  by  allowing 
the  bands  to  slide  across  the  entire  frequency  band  - in 
this  case  the  estimates  would  not  be  independent. 
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Thus,  for  example,  the  first  non-zero  frequency  could 
be  (L  + 2)/2At  and  could  involve  elementary  estimates  2 
through  2L  + 3,  the  next  frequency  could  be  (L  + 3)/2it,  etc, 
Confidence  intervals  for  the  estimates  would  be  obtained 
by  using  the  chi -squared  distribution  and  the  appropriate 
d.f.  \2L  + 1 for  the  center  estimates,  L for  the  end  points). 


{ A - 1 i ) Indirect  Spectral  Estimates 

Suppose  X(j)  = X{jAt),  j=0,  ....  n is  a zero-mean 


process . 


^ 1 n-k-1 

Let  r(k)  = ^ Z X(j)  X(j  + k)  for  k = 0 m (or 

" j = 0 

r(k)  = r(k)  n/(n-k)).  Then 


f(u)  = At  L r(k)e 
k = -m 


■ i 2 ~ u k A t 


, - l/2At  < u < l/2At 


is  another  istimator  of  the  true  spectrum,  based  on  a 
covariance  function  with  maximui..  lag  m.  As  discussed  previously 
the  estimates  must  be  smoothed  with  either  a lag  or  spectral 
window.  In  addition,  the  spectrum  is  only  computed  for 
positive  frequencies  so  the  estimates  are  doubled  to  obtain 
A A m — 1 ^ 

f(u)  = 2At[r(u)  + 2 2 r(k)W(k)  cos(2;TukAt)]  (A. 29) 

k = l 

for  0 < u < l/(2At).  Here  W(k)  = W(kAt)  and  the  ma  X i mum 
lag  is  m. 
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APPENDIX  B 

Cross-Spectral  Analysis 

When  two  or  more  time  series  are  observed,  one  can 
study  their  inter-relationship  by  using  cross-spectral 
analysis.  Subscripts  will  be  used  to  indicate  the  different 
series  and  only  the  case  of  two  series  will  be  discussed. 

Suppose  X^(s)  and  X^Cs)  are  two  zero-mean  time  series 
and  let  X(s)  ==  (X^(s),  X(s)  is  a second-order 

stationary  vector  process  if  the  covariance  matrix,  = 

E(X^(s)  X(s+u))  only  depends  on  u (the  superscript  t indicates 
the  vector  transpose  here).  called  the  cross- 

covariance function.  Note  that 

There  are  two  representation  theorems  that  correspond 
to  Theorems  I and  II  of  Appendix  A,  with  scalers  replaced 
by  vectors  and  matrices. 

In  particular,  for  a second-order  stationary  vector 
process 

X.(t)  = / cos(vt)  oi.(dv)  + / sin(vt)  3-(dv)  (3.1) 

V V J 

where  C0V(aj(dv)  a^(dv))  = C0V{:-j(dv)  5|^{dv))  = Cj^(dv) 

- C0V{:ij(dv)  -^(dv))  - C0V(:..{dv)  ci^(dv))  = qj^.(dv)  and 
all  other  covariances  are  0. 

r(t)  = / e^  ^^  dF(v)  = / cos(t'j)  c(dv)  - / s i n ( tv ) q(dv) 

....  (B.2). 

where  F{v)  is  a matrix  with  Hormitian  non-negative  increments 
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in  V.  c(v)  is  a real  symmetric  matrix  with  entries  Cjj,(v) 
and  is  a real,  skew-symmetri x matric  with  entries 
p^k(v).  c(v)  is  called  the  co-spectral  distribution  and 
q(v)  is  called  the  quadrature  spectral  distribution. 

Pjk  ^ '’jk  C^^(v)] 
is  the  co-herence  spectrum  and 


6jk(^)  “ arctan(-  q^j^  (vO/c^j^  (v))  is  the  phase-spectrum. 
We  will  again  replace  v by  2tu  and  use  capital  letters 
to  denote  the  various  entries,  so  that 


9j^(u)  = arctan[-  Qj|^(u)/C^l^(u)] . 

Once  again  one  can  either  estimate  the  cross  spectral 
densities  directly  by  transforming  the  original  series  or 
by  calculating  cross-covariances  and  then  transforming  the 
cross-covari ances . 


( B . i ) Direct  Estimation 

The  computer  analysis  for  direct  estimation  is  outlined 
below. 

1.  Let  X^(k),k=0,  ....  n-1  and  X^ik),  k=0,  ...,  n-1  be  the 

2 

2 original  series.  Taper  eich  series  with  the  sin  taper. 

2.  Let  Z^(k)  = X.(k-l),  k=1,  ....  n and  fill  out  the  Z-arrays 
to  get  NN  = 2M(2L  +1)  - 1 > n points  if  M + 1 estimates 
are  desired  (each  based  on  2L  + 1 elementary  estimates). 

3.  Transform  each  series  to  get  Z^(k),  i - 1,2,  k = 1,  ...,  NN. 

4.  The  elementary  estimate  of  the  cross  spectrum  is 

F^^d)  = [2  Z^(I)  Z*(I)]At/(.875)NN 


5. 


Average  the  elementary  estimates  over  21  + 1 points 
(except  for  the  first  and  last  estimates  which  are 
obtained  by  averaging  over  L points)  to  get  estimates 
M + 1 estimates,  with  the  middle  estimates  centered 
at  (2L  + 1 )K/2At,  k=l , ....  M. 

6.  The  real  part  of  the  above  average  is  half  of  the 

co-spectrum  and  the  imaginary  part  is  half  of  the 
quadrature  spectrum,  ^ omitting 

the  frequency  argument. 

7.  The  coherence  spectrum  is 

^12  " ^^12  ^12^^’'l  ‘^2 

where  F^  and  F^  are  the  individual  spectra  and  the 
coherence  spectrum  is  calculated  after  averaging  (else 
it  will  be  identically  one). 


The  phase  spectrum  is  5^^  “ arctan  (-  which 

can  be  made  continuous  by  "guessing"  so  that  the  true 
arctangent  rather  than  the  principal  value  is  obtained. 


( B . i i ) Indirect  Spectral  Estimates 

If  is  the  estimated  covariance  with  M lags 

(k=0,  ± 1,  ± 2,  ...,  i M)  then  the  indirect  approach  is 
outlined  below. 

1.  Let  r^{k)  = 

r°(k)  = ■ "i2('k)]. 

2.  Form 

Zl (k)  = r®(k-l ) W(k-1  ) k = l L + 1 

Z2(k)  = r°(k-l)  W(k-1  ) k = l L + 1 

where  W(k)  is  the  lag  window,  (L  £ M). 
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3.  If  the  cross  spectrum  is  desired  :it  frequency  spacing 
1/(N2  At)  from  0 to  l/2At  where  > N,  add  zeroes  to 
the  arrays  to  get  N2  entries.  Once  again,  if  = 2!., 
estimates  are  obtained  at  frequencies  i/2LAt, 

i=0,  1 , , L. 

4.  Fourier  transform,  via  FOURT,  the  2 arrays  to  get 
Z.(lc).  The  smoothed  estimates  of  the  co-  and  quadrature 
spectrum  are 

= 1^2  Re(Z^(k))  - r®(0)]2At 

and  Q^2^^'  " Im(Z2(k))  at  frequencies  ( k-1 )/ ( 2LAt ) , 

k=l , ....  L + 1 . 

( 8 . i i i ) Other  comments  on  Cross-Spectra 

If  both  series  are  subjected  to  the  same  linear  filters 
the  coherence  and  phase  spectra  are  unaffected.  If  a trend 
is  suspected  one  can  difference  the  2 series  and  improve  the 
final  estimates  (see  Jenkins  and  Watts,  Sec.  8.4.5). 

Ore  can  also  re-align  the  2 series  to  get  better  estimates 
of  the  coherence.  The  raw  estimates  of  the  cross-spectrum 
can  be  obtained  without  recalculation  by  multiplying  the 
original  cross-spectrum  by  exp{-  iZ^ru]  where  - is  the  time 
shift.  For  the  estimates  based  on  the  cross-covari ance , new 
even  and  odd  estimates  are  needed  (see  Jenkii.s  and  Watts, 

Sec.  9,3.4). 
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APPENDIX  C 

Variance  and  Covariance  Calculations 

Throughout  this  appendix,  we  assume  we  have  N observa- 
tions, H estimates  and  d degrees  of  freedom.  The  variance 
and ' covari ance  calculations  are  based  on  the  treatment  in 
Hannan,  Chapter  V.  Carats  will  be  used  to  denote  estimates 
and  frequency  will  be  measured  ip  radians.  The  proofs  of 
the  main  theorems  can  be  carried  out  by  using  multivariate 
analysis  (including  complex  multivariate  analysis).  Jenkins 
and  Watts,  Chapter  9,  pp.  372-373  and  appendix  A. 9.1  present 
an  alternative  approach. 

The  following  theorem  is  the  key  theorem  for  practical 
covariance  calculation. 


THEOREM  I 


(Hannan,  Theorem  9,  p.  280),  Under  suitable 


conditions  on  the  moments  of  X(t), 


J!;  s + .P/M)}  = 0, 


(C.l .a) 


for  V ^ 7r/2(mod27T)  . 


for  V ^ 0,  2tt. 


(C.l .b) 


1 im  N 
N- 


g C0V[f.,(v,).fj,,  (v,)]  . 


/or  = Zn . 


89 


Since  most  other  estimates  in  cross-spectral  analysis 
are  functions  of  f.^  and  one  can  derive  the  desired 

covariances  quite  readily. 

For  example,  for  N/H  large  (omitting  the  dependence 
on  v)  we  have  (for  N large) 


(N/M)  VAR(C.j)  = (N/M)  VAR(f..  + f^.) 


» [N/H][VAR(fj  .)  <•  VARIfjj)  + 2 COV  ( f , j , f j 1 


[4K/nd][f,,fjj  (1  - p^.)  . cf./2] 


(C.2) 


Note  that  f 


1 1 


2 is  the  spectrum  of  X^.(t),  and  recall 


that  COV(ZW)  = E(ZW*},  assuming  E(Z)  = E(W)  = 0. 

^ ^ /V  /N  ^ 

Similarly  .since  = (fj^,^  + f^j^)  and 
we  have 

C0V(a,j.£^,)  = C0V[(?.J  * ?..)(?,,  * f,,)] 

' * ^J1  'kk  * ^tk  * ^ek' 

* ^t'ik  'ij  * fjk  'a  * Ut  hj  * 'a 


” f^'ik  *'li  ■ ‘*ik  '^kj  * '■ik  “kj 


On  ‘Iki^^'' 


(C.3) 


VAR(q,.)  = .AR  (f,j  - f,<) 


Id 

^ ^ 

E(f..-2f..  f . . 
' IJ  IJ  IJ 


"P 

f •) 
u ' 


(C.4) 


• A 

'Assume,  without  loss  of  gener'.lity,  that  E(f  ) = 0 for  the 
following  calculations. 
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cov(g,..Q^,)  = - cov  (f,.  - - f,^)) 

' - f|,£  - fij  fkt  * hi 

" ^‘"U  '•kj  ■ '*ik  '^Ij  ■ ‘'u''kj  * '*7k 


COV(Q,.j,C,.)  = 2C,.  Qj^d 


(C.5) 

(C.6) 


COV(C,j,Q,,)  . COV((f,.  * fj,  , ((f^,  - fj,)} 


A «<V  A ^ 


/S  A 


*'kk  - ^k  * <'oi  \i  - *’ji  '^kk* 


f'^kj  ’■ik  * "^ik  '•kj  * “^ki  ‘•ji  * '^jk  '•ki^^*' 


{C.7.a) 


A A 


COV(Q,,.C.j)  .-iCOV«f,,  - f,,)(f,j  . f.,)) 


A A 


A A 


■*='^kk  'tj  - fkk  fij  * fkk  fji  - ^k 


' cov(Cij  g^j) 


(C.7.b) 


The  covariance  of  the  6's  and  p's  can  be  approximated  by 
using  propagation  of  error  - namely  by  expanding  in  a Taylor 
series  and  dropping  higher  order  terms. 

A /S 

Thus,  for  9..  = arctan(-  Q^^/C..)  we  have 

I 0 * V I 0 


«1j  - »'-ctan[6,j/C,j] 


36..  „ 36.. 

3C,j  " 3Q7^  '■‘’ij 

(g,j/(clj  * - ^u' 

(c,j/(c2.  k g?j))(g,j  - Q,j) 


(c.e) 
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. (Qij/(c?j  « var(£,p 

' 2{C..  Q.../(C^  + Q- COV(C  . Q..} 

I J IJ  Ij  IJ  IJ 


Simi  1 an'  ly , 


^ (1  - 


2 x2 


VAR{p‘,)  = (1  - pf-)Vd 


1 J 


(C.9) 


(C.IO) 


Covariances  can  also  be  obtained  by  nul ti plyi ng  the 
appropriate  expansions  and  then  taking  expected  values.  In 
particular  COV(0^ ^ ,0|^^)  will  be  required  for  confidence 
interval  construction  in  the  location  routines. 


C0V(6,.,,^J  ■ [Ou  '’ki  =k:>  ■ '’(j! 


^tj  «k'.> 


* C,j  C,,  COV(Q,.  g,,)]/-:(cfj  ♦ * Qlj-: 

(C.U) 


The  time  lags  estimated  are  linear  combinations  of  the 
0..'s,  if  the  estimates  are  least  squares  estimates.  Thus, 

* J 


N N ^ - 

T,.-  = I e..{p)  p(tv)/  T 
P=1  P=1 


VAR(t.  . ) 


P=1 


p2{iv)2  VAR(T  (p))/(  : p^(Tv)2)2 

» J r>-  1 


C0V(TiiT^^)  = pq  COV(0.j(p).r-^^(q))(j.\)^/[  E^  p^(A.\)2]2 


p=l  q=1 


p = l 

ki 


p = l 


P COV(0..(p),e  {p))(Av)V[  E pMAv)^] 

P-I  ij  Kk. 

N N 

I COV(0  (p).9.  (p))/(  Z p^)^ 


p = l 
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For  spherical  location. 

6^  = (CS), 

6^  = C(S  + '^12' 

5^  ~ C{S  + ^12^ 
are  used. 

Assuming  S and  t..  are  independent,  and 

' M 

of  error,  we  have 

VAR(6^)  = VAR{S) 

VAR(52)  = C^LVAR(S)  + VAR(t^2^^ 

VAR(<52)  = C^[VAR(S)  + VAR{ti3)] 

C0V(6^,52)  = VAR(S) 

C0V(6^.53)  - VAR{S) 

COV{6^.5^)  = C^[VAR(S)  + COV ( v ^ 2 , t ^ 3 ) ] . 

This  covariance  matrix  is  required  in  Append 
For  plane  location,  the  required  covari 


using  propagation 


ix  D. 

a nee  matrix  is 


given  in  Appendix  E. 
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APPENDIX  D 
Spherical  Location 

The  equations  used  for  spherical  location  are  (c.f. 
eqns.  25  - primes  are  dropped  for  typograph i ca 1 convenience). 

(D.1) 


,2  2x2^2 

= r^  + r^  + r3 


x2  / >2  ^ 2 ^ ,2 

62  - (r^  - X2)  + r2  + r3 

^3  " ■ *3^^  ^^*2  ■ -^3^^  '"3 


(D.2) 

(D.3) 


In  these  equations, 

(i) 


5^  = (cS)^,  02  * Ec(S  + -3  " ^ 


12'-’  ’ '3  -13' 

(ii)  c is  the  speed  of  sound, 

( i i i )  S is  the  time  from  initiation  at  the  source  to  its 
reception  at  microphone  K 

(iv)  Tij  (j  = 1,2)  are  time  lags  between  receotion  at 
microphone  1 and  j. 

(v)  , r2,  r3  are  co-ordinates  of  the  source  in  a 

co-ordinate  system  with  origin  at  microphone  1, 
and  x-y  plane  in  the  plane  of  the  three  microphones. 

(vi)  microphone  1 is  at  location  (0,0,0), 

microphone  2 is  at  location  (x2,0,'3)  and 
microphone  3 is  at  location  (x3,y3,0). 

Subtracting  D.2  from  D.l  and  D.3  from  D.l,  we  obtain 

^ (D.4) 


2 2 
6^  - si 


2X2  >"i  - *2 


"^1  ' ^3  ” ^*3  ’"'I  * ^^3  *"2 


’'3  ■ ^3‘ 


(D.5) 
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Hence , 


^1  " ^*1  " ^2  ^ ^2)/(2x2) 

= (6^  ' ^3  ^ *3  "■  >3  ■ ^''3  »'l^/(2y3) 

r3  r^) 


(0.6) 

(0.7) 

(0.8) 


These  locations  have  to  be  converted  back  to  coordinates 
in  the  original  frame  of  reference.  Supoose  the  microphones 
are  at  and  P^,  respectively,  in  the  original  frame 

of  reference. 

(0.9) 

(0.10) 

(D2  • 0^)D.j/(D^  • D.|  ) is  the  projection  of  D2  onto  . 

The  un-norma 1 i zed  co-ordinate  vectors  in  the  new  frame  of 
reference  (the  one  containing  the  microphones  in  the  x-y  plane) 

are 

(0.11) 
(0.12) 
(0.13) 


jiiJWWM 


Microph  t (0,0,0),  microphone  2 is  at  (|D^5,0,0) 

(xj  . = 0)  and  microphone  3 is  at 

((02  • ^ - (Dj  • 

(X3  = (02  • D 3 " < ^2  ^ ^3  = 0. 

Ignorin  slation  by  , the  new  co-ordinate 

vectors  can  ed  as  a linear  combination  of  the  old 

co-ordinate  1*®2’®3 

(0.14) 


(0.74. a) 


The  source  is  at  Z r.U.  and  hence  in  the  o7d  co-ordinate 

.1  = 1 ' ’ 


(0.15) 


(0.75. a) 


where  G is  the  transpose  of  G. 

Element  (i,j)  of  G is  the  j—  element  of  vector  l) . 
so  that 


(D-16) 


(D. 16.a  ) 


'^^nally.  adding  the  translation  p back  in  th 
. , '^1  in,  the  source 

IS  located  at 
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r*  = + 6‘  r. 

next  Obtain  approximate  confidence  limits  for  the 
P-dtcted  position  by  iinearixin,  D.6,  p.7.  and  D.8  with 
respect  to  6^,  s and  6.. 

•■j  ■ r«x  ■■ 

3 

■ J,  “Jk  - E(Jk)h 


(D.17) 


(0.18) 
(0. 18.a} 


j * 

U2.3, 

V/  h e r e 

_ ^''l 
36^ 

<5,/Xj 

^12 

= m 

36^ 

= - 52/x, 

^13 

= . 

553 

= 0 

^21 

. ^^2 
36^  ' 

= - (x^/y^) 

^22 

_ ^^2 
3(52 

/ ^*'1 

^23  ' 

. ^"2 

3(53  ~ 

■ 'V^j) 

®31  " 

, ^^3 
36^  " 

'*1  ■ - ’■7^2, 

^32  ^ 

'362 

*'“1*12  * ''z‘22^^' 

9r, 


(Xj/yjia,, 


(xj/yjJa,, 


®33  ' ■ nj  djj/rj 
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In  vector-matrix  notation 

[r  - E(r)]  = A[6  - E{6)].  (D.lS.b) 

where  A is  a matrix  with  elements  a... 

* w 

The  covariance  matrix  for  r is 
C0V(r,r)  = E{[r  - E(r)][r  - E(r)]^} 

= A C0V(5,5)  A^  (D.19) 

K\ 

The  covariance  matrix  for  r* , the  co-ordinates  of  the 
source  in  the  original  frame  of  reference,  is 

I = C0V(r*  r*)  = G COVCr.r)  G^. 

Hence,  an  approximate  (1  - a)  x 100"  confidence  region 
for  the  source  r*  is  given  by 

(r*  ••  r*}  Z ^ (r*  - r*)  £ X3  {ci/2). 

A simpler  procedure  is  to  use  a Bonferroni  confidence 
region  of  the  form 

■j  ot/  b j 

Since  this  gives  a confidence  box  rather  than  an  ellipsoid. 
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APPENDIX  E 
Plane  Wave  Location 

Here  we  will  use  the  definitions  given  in  Appendix  D. 


Di  = P2 

- Pr  '’2  = 

^3 

^1- 

cos  a^2 

_ ^^12 

■ 

(E.l) 

cos  a^3 

<^^13 

' iPa!! 

(E.2) 

V = 1. 

(E.3) 

V is  the  unit 

velocity 

vector 

of 

the  plane  wave  and  the 

t's  are  time 

lags. 

If  D.  = 

T 

(d-i , d.2. 

d..) 

1 0 

then 

3 

E V.  d 
k = l 

ik  " *^"12 

(E.4) 

3 

E V.  d 
k=l 

2k  " 

(E.5) 

? 2 
k = l 

1 . 

(E.6) 

These  equations  must  be 

solved  to  obtain  (v^,V2,V3) 

- 

''2 

ct^2 

v„ 

a 

<^13 

(E.7) 

Vi  dg^  + 

''2  *^22 

^■^13 

‘^23" 

(E.8) 
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Hence , 


or 


where 


"l  ^ 

[(ct^2  ■ '' 

3di3)d22 

^‘^'^13  ■ 

''3 

^^23 

)d,j 

]/[d„ 

d22 

- c 

’12' 

*21^ 

''2 

[(cTi3  - V 

3‘‘23^'^11 

(cti2 

"3 

^13 

jdj, 

]/[d,, 

d22 

- c 

Il2< 

*21^ 

"1  = 

bii  + b^^'' 

3 

(E. 

,9.{ 

i) 

''2  = 

^2}  ‘^22'' 

3 

(E. 

,9.b) 

*^11 

CcTi2d22 

CTi3d^2^'^t^l  1^22 

- 

'^12 

dj,] 

(E. 

,10, 

.a ) 

^2 

[di2d23  - 

^1 

11^22  ■ 

d^2^21^ 

(E. 

,10 

.b) 

^‘^^13^11  ‘ 

■ ct^2^21^/^‘^11‘^22 

- d 

12^^ 

T 

21  " 

(E. 

,10, 

.c) 

^22 

' f^23^11 

ll‘^22  ■ 

'^1 

2'^21^- 

(E. 

.10, 

.d) 

v2  . v2  - 

^2  3 

1 implies 

(b?2 

+ b^2  ■'■  ^ ) 

''■3  * ''‘’n'’ 

12  * 

21“^ 

22) 

^'3  " 

(b?, 

+ b 

2 

21" 

1) 

= 0 

or 


A V,  + V + C = 0. 
3 3 


Hence , 

V . 


" - 4AC 


B ± 


- AC 


(E.n ) 
2 . ,2 


2A  A 

where  A = b^2  1 » B = ^21^22  ^ ^21  ' 

To  resolve  the  sign  ambiguity,  note  must  be  negative  - 
the  solution  with  the  plus  sign  on  the  radical  appears  because 


the  observed  time  lags  theoretically  could  arise  from  a wave 
traveling  orthogonal  to  the  true  wave.  Hence, 
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Vj  = (-  B - AC)/A 

»l.era  A = bfj  . bfj  * 1.  B - b,,b,j  a b^jb.^  and  C - bf,  * 

3V, 

r j = ■ ^ affr  - A ^i=3/A 


(E-12) 


h ^ 


- AC 


where  3E  _ p, 

3t,-,  " *-^12  ^^^22  ' ^22  cd,,]/(d,,d 


Ln,2  .^22  - =22  Cd,2j/(1„'i22  - djjd^,) 
[2b,,  Cd22  - Bb,,  cdj,  ]/ (d, , d22  - '1,2'<2,) 

■C^  a C2B  jP-  - A '- ]/A 

’■13  3'13  3t,3 


(-  cd,2  b,2  a cd,,  b22)/(d,,d22  - 


— ir  / 

33,3  ' " 21  1 =''l2'>n>'<‘'ll<‘22  - '1l2'<2,)- 

The  covariance  matrix  for  (v,,v„,v„)  is 
r~  12  3 


^2 

^12^22 

^12 

^12^22 

.2 

^22 

^12 

^22 

1 

■ 9v,  , 

- 

OJ 

-*  o 



n 

VAR(t^^) 

+ 1 

and  VAR{v3)  = VAR(r^^)  ^ VAR(r^3) 

Bv-j  3v, 

* 2 C0V(t  T ) ^ 

The  source  is  located  at  - (vpV2.V3)cS  where  S is  the 
time  from  initiation  of  the  source  until  its  reception  at 
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microphone  1.  The  approximate  covariance  matrix  for  the 
source  is 

+ (cS)^  COV(v.v), 

where  COV(v,v)  is  the  covariance  matrix  of  the  's, 

(assuming  S is  independent  of  ^13^‘ 

Once  again  one  can  construct  either  more  exact  ellipsoidal 
confidence  regions  or  box-like  Bonferroni  confidence  regions. 


VAR(T)c‘ 


'^1^2  '*1 


'’2 ''3 


''l''3  ''2''3  ''3 


APPENDIX  F 
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Additional  Remarks  on  Spectra  and  Covariances 

One  can  also  estimate  auto  and  cross  covariances  ty 
transforming  twice  * this  can  actually  save  computer  time 
if  the  covariance  functions  are  required.  (See  Fnochson 
and  Ottnes,  po,  247-248). 

The  procedure  to  obtain  the  auto -covariance  is  as  follows: 

1.  Augment  the  original  series  X(l)  ...,X(n)  with  n zeroes. 

2.  Obtain  the  FFT  of  the  2n-length  sequence,  2(k), 
k=l , . . . ,2n . 

3.  Compute  i-iie  raw  spectrum 

Z(k)  * (At/N)  |Z(k)!^  k*l,  ....  2n. 

4.  Obtain  the  inverse  FFT  of  Z(k). 

5.  Discard  the  last  n points  of  the  inverse  transform 
and  multiply  the  r—  term  fay  (n/n-r)  to  get  the 
covariance  function. 

For  cross-covariances  the  following  procedure  can  be  used: 

1.  Augment  X^(k),  k=l,...,n,  with  n zeroes. 

2.  Obtain  the  FFT  of  both  sequences,  Z^(k)  and  Z2(k). 

3.  Multiply  the  2 sequences  term  by  term,  to  get 
Z^Ck)  = Z^ (k)  Z*(k)  (At/N). 

4.  Compute  the  inverse  transform  of  Z2(k).  Then  the 
first  n points  will  be  C0V(X.jX2)  in  reverse  order 
and  the  last  n points  will  be  C0V(X2Xi)  in  reverse 
order,  after  multiplication  by  (n/n-r). 
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